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SUMMARY 


A mathematical model of an active electromagnetic bearing which includes the actuator, the sensor and 
the control system is developed and implemented in a specialized finite element code for rotordynamic 
analysis. The element formulation and its incorporation in the model of the machine are described in detail. 
A solution procedure, based on a modal approach in which the number of retained modes is controlled by 
the user, is then shown together with other procedures for computing the steady-state response to both static 
and unbalance forces. An example of application shows the numerical results obtained on a model of an 
electric motor suspended on a five active-axis magnetic suspension. The comparison of some of these results 
with the experimental characteristics of the actual system shows the ability of the present model to predict 
its performance. 


INTRODUCTION 


DYNROT is a specialized finite element code for rotordynamic analysis developed at the Department 
of Mechanics of Politecnico di Torino during the last 15 years. Its main features are the use of complex 
co-ordinates for the rotor degrees of freedom linked with flexural behaviour and the subsequent use of mean 
and deviatoric matrices in all cases in which some parts of the machine lack axial symmetry. 

The theoretical bases on which the code is founded have been summarized by the authors in many papers 
(from [1] to [8]). The latest code release (4.1) is based on the only assumption of uncoupling between flexural, 
axial and torsional behaviour of the whole machine; other assumptions as linearity or axial symmetry, which 
were postulated in earlier versions, are no more required. However, the code cannot account for a lack of 
symmetry which interests simultaneously parts of the machine spinning at different speeds: this would make 
it impossible to write the mathematical model in the form of differential equations with constant coefficients 
and would prevent from obtaining exact solutions even in the case of linear systems. 
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The code is implemented in MATLAB (trademark of The MathWorks, Inc.); this allows us to deal with 
complex arithmetic in a straightforward way and improves the portability of the code, including its graphic 
parts, on any hardware (personal computers or mainframes) for which a version of MATLAB exists. 

Aim of the work here described is to develop a "finite element" able to simulate active radial magnetic 
bearings together with their control systems and to implement the solution routines to plot the Campbell 
diagram, the eigen loci and the unbalance response of rotating machinery employing such bearings. This 
allows us to study the effects of the various design parameters of control systems, actuators and sensors on the 
dynamic behaviour of the system. As no assumption on the relative positions of actuators and sensors is done, 
it is possible to study also the dynamic effects of the collocation of such elements. 


MATHEMATICAL MODEL OF THE ROTOR 


Assume that a fixed co-ordinate system xyz is placed with z-axis along the spin axis of the rotor and 
define a set of complex co-ordinates for flexural behaviour as described in [1]. The equation of motion of a 
linear rotating system made by a "rotor" whose spin speed is 0) and a nonrotating "stator", is 


[M] m {q) + ([CL -J® [G]) {q} + ([K ] m - 7 co [C r ]J {q } + 

+[M„]„{«} + 2;0){?})+[C.],|?} +[<:,]/*“'{?} + 

+[K n ] d {? } + ([*,], +;w [C r ]> 2 '“ , {?} = {F C }+ {F n }+<oV r K“' 


where subscripts m and d designate respectively mean and deviatoric matrices; n and r designate nonrotating 

and rotating parts of the system [2] and {q } is the complex conjugate of {q }. 

As the main aim of the present work is the study of the stability "in the small" of the system, nonlinearities 
are not introduced into the model. This is consistent also with the fact that nonlinearities are usually introduced 
in rotordynamics by rolling element or hydrodynamic bearings which are not present in the case of rotors 
running on magnetic bearings. However, if deviations from linearity are present, the system is linearized 
about an equilibrium position. 

The forces exerted by the electromagnetic actuators are introduced in equation (1) through vector {/%.}. 
Such forces are function of the displacement vector {q } owing to the presence of the control system and to 
the characteristic of the actuator itself. They are present also in the study of the free behaviour. 

If the whole system is isotropic with reference to the rotation axis, the free circular whirling can be 
obtained by introducing into the homogeneous equation of motion, which includes the forces due to the 
actuators, a solution of the type 

{q) = {q 0 )e liJ > (2) 


where X is the complex whirl pulsation: its real part is the actual frequency of the whirl motion while its 
imaginary part is the decay rate that must be positive to insure stability. 

The steady-state response to nonrotating forces, as self- weight, can be obtained by introducing a constant 
displacement { q 0 } into equation (1) in which rotating forces {T 7 ,} have been neglected. The steady-state 
response to rotating forces, as unbalance, can be obtained by introducing a constant displacement vector 


320 



{<?} = {?> 


( 3 ) 
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into equation (1) in which only rotating forces {F r } are included together with the forces due to the actuators. 

There is no difficulty to extend the above mentioned solutions to cases in which either the stator or the rotor 
lack axial symmetry, following the lines shown in [2], 

By introducing the mentioned solution vectors { q } into the equation of motion, it can be readily trans- 
formed into an algebraic equation. Note that equation (1) does not include hysteretic damping. Although it 
is possible to introduce all the relevant changes to take it into account, this would compel us to solve the 
eigenproblem linked with free whirling a number of times for each spin speed. In order to avoid this com- 
plication, an "equivalent" viscous damping obtained from the modal analysis of the linear, isotropic, natural 
system will be introduced. Alternatively hysteretic damping can be neglected. This is justified also by the 
reason that hysteretic damping is usually quite small, particularly if compared with the damping introduced 
by the control system, and some approximations on its modelling can be accepted. 


MATHEMATICAL MODEL OF ACTIVE RADIAL ELECTROMAGNETIC BEARINGS 


Nodes and Degrees of Freedom 


A magnetic bearing element of the simplest type must at least include four nodes (Figure 1). Nodes k 
and / are respectively placed at the location of the actuator, the first one on the rotor, the second one on the 
stator. Similarly, nodes p and q are placed at the location of the sensor. If the stator of the machine is not 
introduced into the model, only nodes k and p are included: the actuator is in this case more similar to a 
compliant constraint than to a "spring" element. Note that the co-ordinates of nodes k and l must coincide; 
the same holds for nodes p and q . Each node can be assumed to have only one complex degree of freedom 
related to displacements in direction perpendicular to the spin axis. The complex displacement and the 
corresponding force acting on the node are 


z=x+jy ; F=F x +jF y . 


( 4 ) 
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Figure 1. Sketch of the magnetic bearing 
element with its four nodes. 
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This type of element does not consider the bending torque which can result when the shaft is not aligned 
with the bearing. It can be subdivided into several components: force actuator with power amplifier, position 
transducer and controller [10]. 


Actuator, transducer and controller 


The actuator is physically made by a magnetic circuit equivalent to four electromagnets, placed at angles 
of 90° from each other. A pair of electromagnets is assumed to be located along x -direction, while the other 
is in y-direction. Their characteristics are assumed to be all equal and the shaft is assumed to be centred in 
a reference position. 

The force exerted on the shaft by each electromagnet can be expressed t*y the following function of the 
current / and air gap t 


Constant AT can be either computed using simplified formulae or finite element simulation, or measured 
experimentally. Equation (5), which does not include the saturation of the pole pieces, is strongly nonlinear. 
To linearize its behaviour it is customary to supply the electromagnets with a constant bias current i 0 to which 

a control current i c is superimposed [1], [9]. 

Assuming that the shaft is centred in the reference position with constant clearance c , the force supplied 
by the pair of electromagnets controlling the position along x-axis is 


F=K 


c -x j 


-K 


(i +i Y 

l, o ' v c 

I c +X j 


( 6 ) 


Equation (6) can be linearized for small displacements x and small values of the control current i c as 


dF x dF x i 0 i 0 

F r ,._ =^ri c +^-x=^K- i c + 4 K -x 
OL dx c C 


x lin 


(7) 


The equation linking the stator-rotor interaction force at the actuator location (nodes k and /) with the 
displacements at the same nodes is then 


\ F x, + j F y,\ 







( 8 ) 


The second term at the right hand side describes the behaviour of a spring element with negative stiffness 
and is linked with the unstable behaviour of the uncontrolled magnetic bearing. 
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The actuator is driven by a power amplifier modulated by the controller output signal. In practice, as 
commonly found in magnetic bearing applications, a transconductance amplifier with current feedback is 
used. It provides the control current i c to the inductive load constituted by the pair of electromagnets. 

The resulting transfer function between the controller output signal v u and the control current i c is 


Us) 

v„(j) 


= G a (s) = 


JT„ + 1 


( 9 ) 


where k a is the stationary gain and l/x a is the power amplifier bandwidth. 

The rotor displacements are detected by inductive proximity probes. The transfer function between the 
displacement from the reference position of the rotor and the sensor output signal v x , including the conditioning 
circuitry, is 


vj£) 

x(s) 


= G(s) = 


51 ,+ 1 


( 10 ) 


where k s is the stationary gain and l/% is the sensor bandwidth. 

The transducer dynamics as well as those of the power amplifier should be considered as parasitic 

phenomena and should be placed as high as possible with respect to the desired closed loop system 
bandwidth. 

The controller is a linear filter and can be given with a fixed Proportional-Integral-Derivative (P/D) form 
or with a user defined transfer function. The transfer function of a realistic PID in parallel form is 


PID{s) = k c 


1 + — + 


sTj 


ST; STj+l 


( 11 ) 


where k c is the stationary gain, T t is the reset time, T d is the prediction time and x d is the time constant of the 

causal pole of the derivative term (certainly present due to the parasitic capacitances in the operational 
amplifiers). 

If T i ~^T d > x d , the previous parallel form can be well approximated by the following series form 


K 

PID{s) ~ -f 

* i 


sT, + 1 s T d + 1 
s sx d + 1 


( 12 ) 


It is easily noticed that if the extra-pole 1 / x d may be arbitrarily adjusted then the PID controller can be 

conveniently thought of as the series of a PI and a lead network. In this way it is possible to adjust the 
attenuation of high frequency disturbances. 

When the controller is defined by the user with an arbitrary transfer function the latest has to be expressed 
in the following polynomial form 
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with m <n 


( 13 ) 


C(5)=- 


b m s m + +b t s+b 0 


s n +a n _ 1 s n * + + a x s+a 0 


The user defined transfer function can be used either alone or in series with the PID controller. 


STATE SPACE EQUATIONS IN MODAL FORM 


Modal equations of the mechanical subsystem 


The equation of motion of the mechanical subsystem is equation (1); all matrices which are included in 
it are symmetrical: the complex co-ordinate approach leads to symmetrical gyroscopic and rotating damping 
matrices. It can be transformed in modal form by resorting to the eigenvector matrix [<D] of the undamped, 

natural, axi-symmetrical system: i.e. the matrix of the eigenvectors of matrix [M]~J [K] m . This modal trans- 
formation is always possible and does not introduce approximations in itself but is unable to uncouple the 
equations of motion. Actually it leads to diagonal mean mass and stiffness matrices while damping, gyroscopic 
and all deviatoric matrices are not diagonal and couple the equations of motion. 

A further coupling can come from the fact that when using magnetic bearings the rotor is usually 
unsupported by other means, the only exception being that of hybrid bearing systems. The stiffness matrix 
of the mechanical subsystem is then singular and the first two whirl frequencies are equal to zero. The 
corresponding eigenvectors are usually neither orthogonal with respect to the mass nor orthogonal with 
respect to the stiffness matrix but are linear combinations of vectors which have these properties. The 
consequence is that, ordering the eigenvalues in ascending order, elements in position 1,2 and 2,1 of the 
modal stiffness and mass matrices do not vanish. Although there is no difficulty to substitute the first two 
columns of the matrix of the eigenvectors with their linear combinations having the required orthogonality 
properties, the fact that the first two (rigid body) modes do not uncouple is not considered important and the 
eigenvectors are used as directly obtained. 

The size of the problem can be reduced in two different ways: by reducing the original, non modal 
formulation using Guyan reduction or by performing the modal transformation using a reduced set of 
eigenvectors. Both reduction schemes yield approximated results that can be quite close to those of the original 
problem if the reduction is performed carefully. 

The strategy here considered is to start with Guyan reduction to reduce the size of the eigenproblem 
needed to perform the modal transformation. Since Guyan reduction is standard in DYNROT code [1], it 
will not be dealt with here any further. 

The use of a reduced set of eigenvectors can give way to the usual spillover problems and must be 
performed with care. Here the modes are not only coupled by the control system, as usual, but also by the 
characteristics of the mechanical elements and, in particular, by the gyroscopic matrix and, if the system has 
no axial symmetry, by the deviatoric matrices. Several tests performed using a different number of eigen- 
vectors did show that very accurate results can be obtained by using a number of modes slightly higher than 
the number of critical speeds of the free -rotor system encountered within the working range of the machine. 

When in doubt, the user can perform the computation retaining a different number of modes and compare 
the results obtained. 
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Isotropic machine 


The state equations are here obtained in explicit form only for the case of an isotropic machine; their 
extension to the more general case is straightforward. 

By neglecting the deviatoric matrices and introducing explicitly the forces {F c } exerted by the actuators, 
equation (1) yields the following modal equation 

m {Tl} + ac] -yco [G]) (li) + ([K] -j ( o [C r ]) {ri} = 

= {F c } + {F n } + <o 2 {Fr}e j(a ' , (14) 


Vector {r|} is the modal complex co-ordinates vector ((T)} = [<1>] {< 7 }); its real part refers to the 
xz -plane, while its imaginary part to the yz -plane. 

The mass and stiffness matrices ([M] = [d>] r [M] [d>] and [Af] = [<J>] r [Af] [O]) are diagonal and coincide 
with the conventional modal matrices linked with the vibration of the system in the xz -plane. 

The damping matrices (global damping [C] = [d>] T [C] [<I>] and rotating damping [C r ] = [<I>] r [C r ] [O]) are 
not diagonal except for the case of generalized proportional damping. However, as usually the mechanical 
subsystem is lightly damped, they can be approximated by simply neglecting their elements outside the main 
diagonal. 

The gyroscopic matrix ([G ] = [d>] r [G ] [<!>]) is not diagonal and, except for the case of very slender rotors 
which exhibit a very small gyroscopic effect, cannot be approximated by simply neglecting the elements 
outside the main diagonal. 

Vectors {F n } and {F r } are respectively the complex modal force vectors linked with stationary and 

unbalance forces. Their real parts refer to the xz-plane (stationary reference system or rotor-fixed reference 
respectively), while their imaginary parts refer to the yz -plane. 

Vector (F c } is the complex modal control force vector ({A 7 ,. } = [<J>] r [S c ]{/\.}). Vector {F c }contains the con- 
trol forces in complex modal form: the real part refers to the xz -plane, the imaginary part refers to the yz -plane; 
the number of its elements is equal to the number of actuators. Matrix [SJ has a number of lines equal to the 

number of degrees of freedom of the system and a number of columns equal to the number of nodes related 
to the actuators. All its elements are zero except for the elements ih the column related to each bearing and 
the line related to the degree of freedom of the rotor at the actuator location (displacement at the A:-th or /- th 
node, see Figure 1) whose value is 1. 

Expressing the force vector {F c } as a function of the control currents and displacements through equation 
( 8 ), the modal control force vector {F f } can be written as 

(F.) = -[4>f [SJ [A" ] {!,}- [4-f [S t ] IK.} [<i>] {til = -[£,] (U - \K.} {Tl} , (15) 


where [AT,] is a matrix with a number of rows equal to the number of nodes related to the actuators and a 
number of columns equal to the number of the latter; it is simply obtained from the first term at right hand 
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side of equation ( 8 ). In a similar way, [K u ] is a square matrix with the same number of rows as [A',] obtained 
from the second term at the right hand side of equation ( 8 ); it is nonpositive defined and expresses the open 
loop negative stiffness of the actuators. 

Equation (14) can thus be expressed in the form 

{^} + [M ]" 1 «C] -yco [G]) {ii) + [M]' (M + [KJ ~j<» [C r ]) {t\} = 

= -[Mf [*,] (U + ffif {FJ + wWr 1 i F r }e im ‘ . 

The complex output vector {yj containing the displacements at the sensor locations can be expressed 
as a function of the modal co-ordinates as 

{y,}=[s.H<t>i(m . < n > 

where matrix [SJ has a number of lines equal to the number of sensors and a number of columns equal to 

the number of degrees of freedom of the system. All its elements are zero except for the elements in the row 
related to each sensor and the column related to the degree of freedom of the rotor at the sensor location 
(displacement at the p-th node, see Figure 1) whose value is 1 and to that in the same row and the column 
related to the degree of freedom of the stator at the sensor location, if any (displacement at the 17 -th node), 
whose value is - 1 . 


Open loop state space equation 


By introducing a state vector including the modal velocities and displacements, the open loop state 
equation of the mechanical sub-system is 


| (z } = [*] (z } + [B c ] (O + [B n ] {F n } + a > 2 [B r ] { F 

1 W = [CHz} + [®HU . 

where 


(z} = 

[B c ] = 


{m 

.{nl 


-[Mf 1 [AT,] 


[*] = 


-[Mf 1 ([C] -y to [G]) - [Mf ([K] + [K u ] -j (0 [C r ]) 


[/] 


[ 0 ] 


, [B n ] = [B r ] = 


[M] 

[ 0 ] 


-1 


[ 0 ] 


, [C}= [0] [5J[<D] , 


[®] = [ 0 ] 


(18) 


(19) 


The manipulable command signal {/ c } and the non manipulable signals (Fn), co 2 {F r }e ja * can be organized 
in a unique term. In this way equation (18) can be rewritten in the following more compact form 
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{z} = [A]{z} + ['B]{u} 

•W = [C]{Z}+[©]{U , 


(20) 


where 


{«} 


' (U ' 


{F n } 

> , [«] = 

0> 2 {F r }e M 



[B J 


[BJ W r ] 


( 21 ) 


Note that the dynamic matrix of the system and all vectors are complex. An equivalent formulation 
employing real co-ordinates is possible, but it would involve matrices of double size. 


Closed loop state space equation 


While the (linearized) force actuator is already taken into account in the open loop state equation of the 
isotropic machine, the controller, the power amplifier and the position transducer are usually described as 
single-input single-output transfer functions in the frequency domain. Thus the Laplace domain descriptions 
given in the previous section must be converted to a time domain form to be compatible with the rotor model 
equation (1) and open loop state equations (20). The desired form is the state space formulation obtained by 
a minimal realization (of equations from (9) to (13)) in the controllable canonical form. 

Since the open loop state variable is in complex form also, the single-input single-output realization for the 
controller, the power amplifier and the position transducer should be thought of accordingly. It is to be noted 
that the controller may be directly designed as a multiple-input multiple-output system (i.e. by pole placement 
techniques) and represented in complex matrix form. 

Given the state equation (20) of the plant and given the state space realization of controller, actuator and 
transducer with the following quadruples 
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the closed loop equations are obtained through standard series and then feedback connections as shown in 
Figure 2. 


|F n (t)i c^F r ie jut 



Transducer 


Figure 2. Closed loop blocks and 
signals diagram. 
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Here (v r (f)} is the position reference signal, {v c (0) the command signal to the power amplifier, {/ c (0) 
the control current (then force), {F,,}, co 2 {/ 7 r }^ /OM the rotor stationary and unbalance modal forces, { } 

the displacements at sensor location and { v,(r ) } the sensors output. 

In the regulation problem the reference signal may be put to zero and the closed loop system^output 
{^(/)| expressed in function of the rotor stationary force {F n ) and rotor unbalance modal force co 2 {/\}e'"' 


{Cl + 

{.v,} = [C CL 1 {Q + [Vcl 1 


(23) 


where {£} is an extended state vector containing all modal states (z } of the system plus the realized state of 
the control loop. 


Campbell diagram and eigen loci 


The Campbell diagram and the velocity eigen locus are two ways of representing in graphic form the 
variation of the whirl frequencies and decay rates when the spin speed changes. They are essentially equivalent 
but, as the first is more common with rotordynamics specialists while the second with control scientists, 
the DYNROT code is designed to supply both. 

Usually when computing the Campbell diagram the solution of the equation of motion is assumed to be 
of the type (2) and consequently the whirl frequency is the real part of X while the decay rate is its imaginary 
part. The real and imaginary parts are plotted separately as functions of the spin speed 0). 

When plotting the velocity eigen locus a solution of the type 


{q} = {q 0 }e s ‘ 


(24) 


is assumed: the whirl frequency is thus the imaginary part of s while the decay rate is its real part changed 
of sign. The imaginary part is then plotted as a function of the real part in the complex plane (0) does not 
appear explicitly). 

As the system is assumed to be linear, the Campbell diagram depends neither on unbalance nor on static 
external forces. By introducing the solution (24) into the closed loop state equation in modal form (23), in 
which the external forcing functions have been neglected, it is readily transformed into the algebraic equation 

( \ < 25 > 


328 


Matrix [\ L ] is in this case complex and, as a consequence, the eigenvalues are not conjugate. This 

feature, which is uncommon in the solution of the eigenproblem related to state equations, is directly linked 
with the use of complex co-ordinates and with their physical meaning: the real and imaginary parts are the 
generalized displacements (or modal displacements) in xz- and yz - planes and the sign of the imaginary part 
of s has a physical meaning, indicating whether the circular whirling occurs either in forward or in backward 
direction. As a consequence of the gyroscopic effect the whirl speeds for forward and backward whirling are 
different and the solutions for s cannot be conjugate. A different situation occurs in the case of unsymmetrical 
systems: elliptical whirling is here obtained as the combination of a forward and a backward circular whirling 
at the same frequency, and then solutions with an equal imaginary part but with anoppositesign must be present. 

Apart from plotting the Campbell diagram and the decay rate plot as a function of the speed and the 
velocity eigen locus at varying speeds, a routine which computes the closed loop eigen locus at fixed speed 
with varying the overall gain of the feedback loop has been developed. 


Steady-state response 


The response to unbalance and to constant nonrotating forces can be computed as 


M/l-Mcd {Q = co 2 [fi r ]{F r } 




(26) 


For each value of the speed spin to, equation (26) can be readily solved in {£<,}. After applying the inverse 
modal transformation the same solutions yield the response of the system. Note that in general and also in 
the case of response to static forces the dynamic matrix of the system depends on the spin speed. 


NUMERICAL EXAMPLE 


The five-active axes magnetic suspension described in [11] was used as a test case for the present 
mathematical model and then for the code used to implement it. 

The rotor has been modelled with 22 beam elements and 3 concentrated mass elements. Two magnetic 
bearing elements have been added and the stator was considered as' a rigid body. Consequently the nodes of 
the type / and q of Figure 1 were directly grounded. The model has a total number of 46 degrees of freedom, 
reduced to 23 by using the Guyan algorithm considering as master degrees of freedom all complex 
translations. The damping of the material constituting the shaft has been neglected in the model. 

A first analysis of the open loop system was performed. The first 4 modes corresponding to the critical 
speeds are reported in Figure 3b together with a sketch of the model of the shaft (Figure 3a). The first two 
modes are rigid-body modes and the other two correspond to critical speeds at 25780 and 54890 rpm. 

The computation was repeated using a PID controller on each axis; the results for the first 6 modes are 
reported in Figure 3c. Here the first two modes are rigid-body modes at 2990 and 3880 rpm. They are well 
damped, the decay rates being 17.4 and 28.9 1/s, respectively. Two additional modes linked with the control 
system can then be found at 7290 and 7350 rpm; they are however very much damped, with a decay rate of 
20930 1/s and cannot be expected to be seen experimentally. The first two deformation modes of the rotor 
occur at 25860 and 54810 rpm, essentially unaffected by the presence of the controller. 
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Rotor with five-active axes magnetic suspension 



Figure 3. a) Sketch of the FEM model of the rotor, b) 
First 4 modes corresponding to the critical speeds of the 
open loop system, c) First 6 modes corresponding to the 
critical speeds of the closed loop system with PID con- 
troller. (Screen dumps produced by DYNROT code). 
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Figure 4. Campbell diagram: whirl speed and decay rate as functions of the spin speed. (Screen dumps produced by 
DYNROT code). 
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Figure 5. Velocity eigen loci at varying spin speed 
in the range 0 - 5000 rad/s. (Screen dump produced 
by DYNROT code). 


Root locus 



The decay rate corresponding to the last two modes is still positive but quite small (1.86 and 1.34 1/s). 
Other critical speeds, well outside the working range, have negative decay rates showing that at those speeds 
the system would be unstable. 

The Campbell diagram is shown in Figure 4. It has been obtained by considering only 5 modes as adding 
further modes would not change the results in the field of speeds and frequencies of interest. Note that the 
whirl frequencies are not much influenced by the spin speed, as was easily predictable from the geometrical 
configuration of the rotor. Actually in the case of rotors with small gyroscopic effects as the present one the 
whirl frequencies can be assumed to be independent from to and the information gained from the Campbell 
diagram is not as important. A greater dependence from the spin speed is shown by the decay rate. The 
controller modes mentioned above are reported in the plot of the whirl frequencies but not in that of the decay 
rate, as they are out of scale. 

The velocity eigen loci at varying spin speeds are shown in Figure 5. The plot has been obtained for speeds 
ranging from 0 to 5000 rad/s, with increments of 1000 rad/s: a finer pitch is considered useless as the results 
are not much affected by the spin speed. Note that in the present case the velocity eigen loci are not symmetrical 
with respect to the real axis: in the present case, negative values of the imaginary part of s designate backward 
whirling and positive values of s indicate forward whirling. They are always different even if in the present 
case, owing to the low gyroscopic effect, this difference is not great. 

The unbalance response has then been computed. An eccentricity of 1 pm of rotor of the electric motor 
has been assumed, which corresponds to a quality balancing grade G 3.6 at 35000 rpm, following the 
ISO 1925 standard. 

The amplitude of the orbit at the location of the motor (node 5) and of the sensor of the radial bearing at 
the right in Figure 3a (node 16) are shown in Figure 6. The amplitude at the rigid-body critical speeds is far 
smaller than that at the first deformation mode; the two controller modes are completely damped out. At the 
motor location the two rigid-body critical speeds are well separated and an ample zone in which the amplitude 
is very small is visible between the rigid-body and the deformation critical speeds. Self-centring occurs in 
most of the working range above the rigid-body modes, except near the third critical speed. The response at 
the sensor location shows slightly different patterns: the rigid body critical speeds are not well separated, 
tending to merge into a single peak. 
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Figure 6. Unbalance response. Amplitude of the orbits (in m) at the centre of the motor (node 7) and at the sensor of 
the bearing at the right (node 16) caused by an eccentricity of 1 pm of the rotor of the motor (balancing quality grade 
G 3.6 at 35000 rpm). (Screen dumps produced by DYNROT code). 

The natural frequencies of the system at standstill obtained from the Campbell diagram of Figure 4 are 
compared in Table I with the experimental values obtained by exciting the system with an instrumented 
hammer at the centre of the shaft and measuring the acceleration in the same point. The values obtained from 
the Campbell diagram in correspondence with the frequencies due to the controller have not been reported. 

The accordance is very good, particularly considering that the comparison has not been made between 
exactly correspondent quantities. The experimental values correspond to the peaks of the impulse response 
while the numerical ones are the damped eigenvalues of the system. While in the case of lightly damped 
modes no noticeable difference exists between them; in the case of modeswith higher damping, the comparison 
can be at best indicative. The first two natural frequencies merge into a single value; a similar effect can be 
seen in the numerical results of Figure 6. 

Table I. Experimental Values of the Lowest Natural Frequencies (Hz) Compared with the Values Obtained from the 
Campbell Diagram of Figure 4 at Standstill 


Experimental 

57 

420 

884 

Numerical 

47.9 

60.8 

419 

855 


CONCLUSIONS 


A mathematical model for the dynamic simulation of rotating machines with active electromagnetic 
bearings based on the FEM and the use of complex co-ordinates has been described in detail. The main aim 
of the model and of the code based on it is the computation of the Campbell diagram to study the stability 
in the small and of the steady-state responses to unbalance and static forces. 
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At present the model is limited to machines possessing axial symmetry; however, this limitation is not 
intrinsic to the model used and can be easily avoided by simple changes to the solution algorithms. Another 
limitation which can be easily eliminated is the assumption that the controller operates on the two axes of 
each bearing and on the various bearings as a number of separated SISO units. There is no difficulty in 
introducing a complex transfer function to allow coupling between the control in xz- and yz- planes or to 
couple together the various controllers. 

An assumption which on the contrary is actually needed is that of linearity. It can be dropped by introducing 
nonlinear elements but this would make it impossible to obtain Campbell diagrams or eigen loci and the only 
possible computation would then be that of the unbalance response. However, the small values of the 
amplitudes which can be obtained using magnetic bearings guarantee that the system remains in most cases 
within limits in which the linearity of the model does not introduce unacceptable approximations. 

The possibility of introducing an accurate modelization of complex rotors allows us to obtain accurate 
results for practical cases and in particular allows us to study the effects of the relative position of the 
sensors and the actuators (collocation problem). 

An example has shown the types of results which can be obtained from the routines developed from the 
equations here described and incorporated in the DYNROT finite element code. The comparison with some 
experimental results shows that a very good accuracy can be obtained. 
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